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EXECUTIVESUMMARY 

The  msin  objective  of  the  program  (F49620-93- 1-0513)  was  to  determine  if  and  where 
coherent  transient  technology  can  be  used  in  future  memory  and  signal  processing  devices,  to 
develop  techniques  for  implementing  practical  devices,  and  to  demonstrate  the  feasibility  of  the 
applicable  devices.  The  approach  encompassed  six  main  areas:  1)  analysis  and  simulation, 

2)  evaluation,  3)  experimental  validation  of  theory,  4)  material  exploration  and 
development,  5)  feasibility  demonstrations,  and  6)  collaborative  interactions  with  material 
growers  and  system  designers. 

Almost  all  the  program’s  findings  have  been  or  will  be  published  in  archival  journals  or 
proceedings.  The  program  resulted  in  9  published  and  4  to  be  submitted  journal  papers,  7 
refereed  conference  publications,  5  invited  talks,  several  other  conference  and  workshop 
presentations,  and  2  patents.  For  those  findings  describe  elsewhere,  brief  summaries  of  the 
work  accomplished  are  given  along  with  references  to  the  relevant  published  articles  or 
proceedings.  The  work  accomplished  under  this  program  was  augmented  by  an  ASSERT  grant 
F49620-95- 1 -0468  and  DURIP  equipment  F49620-95- 1  -0477. 

The  significant  accomplishments  of  the  program  include: 

1 )  a  theoretical  evaluation  of  the  performance  of  the  coherent  transient  continuous  optical 

processor. 

2)  the  demonstration  of  the  continuous  correlator  showing  that  correlation  and  convolution  of 

signals  can  be  done  well  beyond  the  homogeneous  decay  time  limit. 

3)  an  analysis  of  the  performance  of  coherent  transient  memory  systems  based  on  physical 

limits. 

4)  the  development  of  a  novel  concept  in  optical  routing  and  processing  that  combines  the 

temporal  processing  capabilities  of  coherent  transients  with  the  spatial  processing  of 
volume  holography. 

5)  the  proposal  and  demonstration  of  novel  techniques  for  header/data  isolation  based  on 

spatial-spectral  holography. 

6)  the  analytical  development  of  scaling  relations  that  relate  material  properties  to  the  system 

data  bandwidth. 

7)  the  demonstration  of  the  replacement  of  the  single  brief  reference  pulse  used  in  coherent 

transient  processors  with  two  long,  frequency  chirped  pulses. 

8)  the  development  and  demonstration  of  techniques  to  compensate  for  the  deleterious  effects 

of  homogeneous  dephasing  in  coherent  transient  optical  memories  and  processors. 

9)  The  development  and  demonstration  of  the  application  of  spatial-spectral  holography  to 

true-time  delay  processing  for  phased  array  radar  transmitters  and  receivers. 

10)  An  analysis  of  the  enhanced  efficiencies  that  are  achievable  in  inverted  media  and  in 
optically  thick  media. 

11)  The  demonstration  of  the  use  of  spatial  spectral  holography  for  frequency  division 
multiplexed  communication,  interconnection,  and  routing. 

12)  Further  development  of  an  optical  coherent  transient  signal  simulator. 
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Performance  Evaluation  of  the  Coherent  Transient  Continuous  Optical  Correlator 

In  collaboration  with  John  A.  Bell  at  Boeing,  a  theoretical  evaluation  of  the  performance  of 
the  coherent  transient  continuous  optical  processor  was  completed.  The  analysis  considered 
the  effects  of  coherent  saturation,  incoherent  saturation,  gating  efficiency,  pattern  pulse 
compression,  homogeneous  decay,  cooperative  emission,  spatial  beam  overlap,  spatial 
crosstalk  between  adjacent  beams,  diffraction,  and  local  thermal  heating.  The  main  inputs  to 
the  analysis  are  the  wavelength,  the  inhomogeneous  linewidth.  the  homogeneous  linewidth,  the 
number  density  of  absorbing  centers,  gating  efficiency,  index  of  refraction,  and  detector 
efficiency.  Additional  inputs  involve  the  spectral  and  spatial  characteristics  of  the  pattern  and 
data  streams  and  the  acceptable  level  of  nonlinearities.  Given  the  number  of  detected  photons 
per  bit  needed  to  achieve  a  desired  signal  to  noise,  the  analysis  yields  the  required  spot  size  and 
maximum  pattern  density  for  a  given  set  of  input  parameters.  For  example,  for  a  material  with 
reasonably  obtainable  characteristics  the  correlator’s  expected  data  rate,  time-bandwidth 
product,  and  pattern  storage-density  are  respectively  5  gigahertz,  16000,  and  1 10,000  patterns 
per  square  centimeter.  Higher  performance  depends  on  developing  higher  performance 
materials.  The  calculation  also  yields  the  optimal  oscillator  strength  and  path  length  through  the 

crystal.  The  work  was  published  in  Applied  Optics  33,  1538-1548  (1994)  and  resulted  in  U 
S.  Patent  No.  5,239,548. 

Coherent  Transient  Continuous  Optical  Correlator  Demonstration 

In  collaboration  with  Miao  Zhu  and  C.  Michael  Jefferson  at  IBM,  it  was  demonstrated  that 
correlahon  and  convolution  can  be  done  with  coherent  transients  well  beyond  the  homogeneous 
decay  time  limit.  The  basic  concept  that  was  demonstrated  was  that  once  the  spectral 
interference  grating  is  stored  in  the  ground  state,  it  can  be  continuously  interrogated  provided 
the  transirion  is  not  saturated.  This  allows  the  material  to  act  as  a  passive  spatial-spectral  filter 
with  application  in  signal  processing,  database  memory,  and  routing  applications.  In  the 
demonstration,  a  3 120-bit  long  phase  encoded  data  stream  was  correlated  with  a  13-bit  phase 
encoded  barker  code.  Even  though  the  data  stream’s  duration  exceeded  the  upper  state  lifetime 
of  the  excited  transition  (Thulium  doped  Y ttrium  Aluminum  Oxide),  the  output  signal  intensity 
maintained  level  and  the  output’s  fidelity  was  excellent.  This  demonstrates  ability  to  perform 
continuous  real-time  processing  of  phase  and  amplitude  encoded  optical  data  streams  using 
coherent  transient  techniques.  The  results  were  presented  as  an  invited  talk  at  Spectral  Hole- 
Burning  and  Related  Spectroscopies:  Science  and  Applications,  August  24-25  1994  and 
published  in  Optics  Letters  2  0,  2514-2516  (1995). 

Physical  Limits  Of  Coherent  Transient  Memories  With  Two-Dimensionall  Addressing 
In  collaboration  with  Thomas  W.  Mossberg  of  the  University  of  Oregon,  an  analysis  was 
made  of  the  performance  of  coherent  transient  memory  systems  in  which  multiple  bits  are 
temporally  stored  at  each  spatially  addressed  laser  spot.  This  approach  differed  from  the 
approach  used  in  an  earlier  papeifW.  R.  Babbitt,  Proc.  SPIE  2026, 532-542  (1993)),  which  is 
similar  to  the  approach  used  in  the  evaluation  of  the  continuous  correlator.  In  the  earlier 
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memory  and  continuous  correlator  analyses,  a  particular  implementation  and  set  of  performance 
criteria  were  assumed,  which  yielded  an  overall  efficiency  of  the  process  that  was  particular  to 
the  chosen  inputs.  In  this  paper,  no  particular  implementation  was  assumed.  Instead,  the 
overall  efficiency  of  the  process  is  an  input  parameter.  Based  on  consideration  of  the  physical 
limits  involved  in  the  two-dimensional  storage  process,  limits  were  placed  on  the  maximum 
obtainable  areal  storage  densities  for  coherent  transient  memories.  For  a  material  that  allows 
one  million  bits  per  laser  spot,  we  fold  a  20,000-fold  increase  in  areal  storage  over  traditional 
2D  optical  memories.  These  are  the  physical  limitations  and  known  practical  considerations 
and  implementations  may  not  enable  these  densities  to  be  obtain.  However,  novel 
implementations  are  continually  being  developed  that  push  the  technology  closer  to  the  physical 
limits.  In  addition,  higher  performance  can  be  obtained  using  volume  holographic  storage 
capabilities  of  coherent  transients.  The  physical  limits  of  coherent  transient  volume  storage 
have  yet  to  be  analyzed,  but  initial  calculations  indicate  significantly  higher  areal  densities  can 
be  achieved  with  relaxed  absorber  concentration  requirements.  The  results  of  the  physical 
limits  of  quasi  two-dimensional  storage  were  published  in  JOSA  B  1 1,  1948-1953  (1994). 
Spatial-Spectral  Holographic  Routing  Of  Optical  Beams 

In  collaboration  with  Thomas  W.  Mossberg  of  the  University  of  Oregon,  a  novel  concept 
in  optical  routing  and  processing  that  combines  the  temporal  processing  capabilities  of  coherent 
transients  with  the  spatial  processing  of  volume  holography.  This  new  concept,  referred  to  as 
spatial-spectral  holography,  has  applications  in  code-division  multiplexed  routing,  packet 
header  recognition,  header/data  isolation,  and  contention  resolution.  A  patent  application  has 
been  filed  jointly.  The  work  resulted  in  several  conference  presentations  and  in  the  paper 

“Spatial  routing  of  optical  beams  through  time-domain  spatial-spectral  filtering”  (Optics  Letters 
20,910-912). 

Optical  Coherent  Transient  Header/Data  Isolation 

Novel  techniques  for  header/data  isolation  based  on  spatial-spectral  holography  were 
proposed  and  demonstrated.  Using  these  techniques,  ultra-high  bandwidth  optical  information 
stream  can  be  parsed  into  two  segments.  Either  segment  can  be  recalled  separately  (isolated 
from  the  other).  As  an  example,  an  information  packet  with  a  header  and  a  data  segment  can  be 
parsed.  The  isolated  header  can  be  processed  and  then  the  stored  data  segment  can  be  routed  or 
processed.  The  novel  header  and  data  isolation  techniques  and  their  experimental 
demonstrations  were  presented  at  HBRS’%and  in  Optics  Letters  21, 71  (1996). 

Scaling  Of  Material  Properties  With  Data  Bandwidth 

One  important  finding  in  the  evaluations  of  coherent  transient  systems  is  that  as  the 
bandwidth  of  devices  are  pushed  to  higher  bandwidths,  the  material  properties  that  relate  to 
rates  must  also  be  proportionally  increased.  Most  significantly,  the  optimal  oscillator  strength 
scales  with  the  data  bandwidth.  This  may  eliminate  rare  earth  f-f  transition  from  consideration 
in  terahertz  devices.  This  finding  assumes  a  constant  number  of  detected  photons  per  bit  to 
achieve  a  desired  bit  error  rate.  As  bandwidths  increase,  a  higher  number  of  photons  per  bit 
may  be  required.  Under  these  conditions,  the  bandwidth  and  material  rates  won’t  scale 
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proportionately,  but  the  trend  will  still  be  towards  higher  oscillator  strengths.  These  findings 
have  been  presented  and  discussed  at  the  several  workshops  on  applications  of  spectral 
holebuming. 

Chirped  Reference  Pulses  In  Coherent  Transient  Processing 

The  replacement  of  the  single  brief  reference  pulse  used  in  coherent  transient  processors 
with  two  long,  frequency  chirped  pulses  was  demonstrated.  Himinating  the  intense  brief  pulse 
requirement  makes  programming  coherent  transient  processors  considerably  more  practical. 

As  proposed  in  the  original  patent  on  coherent  transient  processing,  the  first  brief  pulse  in  triple 
product  correlation  sequence  could  be  replaced  by  two  frequency  chirped  pulses.  We  have 
both  analytically  and  experimentally  shown  this.  We  have  also  performed  calculations  that 
show  that  another  previously  suggested  implementation  will  not  work,  but  that  three  novel 
implementation  we  have  proposed  we  work.  One  of  these  novel  implementations  was 
demonstrated  as  well  as  the  other  previously  proposed  implementation.  The  demonstration 
results  were  presented  at  OSA  Annual  Meeting  1994,  Photonics  in  Computing  II,  and  in 
Applied  Optics  35, 278  (1996). 

Compensation  for  homogeneous  dephasing  in  coherent  transient  systems 

Techniques  to  compensate  for  the  deleterious  effects  of  homogeneous  dephasing  in 
coherent  transient  optical  memories  and  processors  were  developed  and  demonstrated.  By 
tailonng  the  pattern  pulses  that  program  the  coherent  transient  material,  distortions  of  output 
signals  can  be  significantly  reduced,  allowing  time-bandwidth  product  or  storage  capacity  of 
coherent  transient  devices  to  be  increased  by  up  to  an  order  of  magnitude.  These  techniques 
may  aid  in  the  practical  implementation  and  enhanced  performance  of  coherent  transient 
routers,  processors,  and  memories.  Trade-offs  with  signal  size  must  be  considered  in 
evaluating  the  overall  performance  of  a  system  implementing  homogeneous  decay 
compensation.  The  techniques  are  described  and  demonstrated  in  the  proceedings  of  CLEO’96 
and  HBRS’96  and  in  Optics  Communications  128,  136  ( 1996). 

Coherent  transient  true-time  delay  regenerators  and  processors 

The  application  of  spatial-spectral  holography  to  true-time  delay  processing  for  phased 
array  radar  transmitters  and  receivers  has  been  developed.  Coherent  transient  true-time  delays 
have  the  ability  to  simultaneously,  and  continuously  process  several  hundred  delays  over  a 
broad  bandwidth  with  fine  temporal  resolution.  This  will  enable  wide  angle  (>  45  degrees) 
beamsteering  of  wide-bandwidth  (>  lOGHz)  for  multi-element  (>  1000)  arrays. 

Multicasting  and  simultaneous  delay  and  temporal  processing  capabilities  are  additional 
features.  The  concepts  and  their  demonstration  have  been  presented  at  several  conferences, 
including  HBRS’96  and  in  Optics  Letters  21, 1 102  (1996).  Various  techniques  for 
programming  the  coherent  transient  true-time  delay  regenerator  using  frequency  chirped  pulses 
are  described  in  a  manuscript  being  submitted  to  Optics  Letters.  The  role  of  non-linearities  in 
the  material  excitation  of  the  delay  timing  is  described  in  a  manuscript  being  prepared  for 
submission  to  Physics  Review  Letters. 
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Echoes  in  an  Inverted  Media  and  Optically  Thick  Media 

The  enhanced  efficiencies  that  are  achievable  in  inverted  media  has  been  shown  via 
computer  simulation.  The  most  significant  result  is  that  data  recall  efficiencies  much  greater 
than  unity  are  achievable  in  optically  thick  inverted  media  while  maintaining  high  fidelity. 
Enhanced  efficiencies  are  critical  to  making  coherent  transient  devices  practicable.  The  results 
for  the  inverted  media  were  described  at  the  HBRS’96  conference  and  in  a  paper  published  in 
Mol.  Cryst.  Liq.  Cryst.  29 1, 269  (1996).  A  comprehensive  study  of  the  efficiencies  of  data 
pulse  recall  in  non-inverted  optically  thick  media  is  being  continued  under  our  new  AFOSR 
grant  and  a  paper  is  being  subimtted  to  Phys.  Rev.  A.  The  most  sigmficant  result  here  is  that 
efficiencies  approaching  unity  are  achievable  in  even  non-inverted  media,  which  is  counter  to 
the  conventional  wisdom.  Experiments  to  verify  these  results  are  in  progress. 

Spatial  Spectral  Holography  Frequency  Division  Multiplexing 

The  proposed  use  of  spatial  spectral  holography  (SSH)  for  frequency  division 
multiplexing  (FDM)  was  demonstrated.  In  contrast  to  volume  holographic  EDM,  where  the 
efficiencies  are  inversely  proportional  the  square  of  the  number  of  channels,  the  efficiency  of 
SSH  FDM  routing  is  independent  of  the  number  of  channels.  For  a  1000  channel  FDM 
system,  the  10-25%  channel  efficiency  of  SSH  FDM  is  extraordinary  compared  to  efficiencies 
of  on  the  order  of  ( 1/10(X))2  per  channel  achievable  with  volume  holography.  With  channel 
bandwidths  over  a  gigahertz,  SSH  FDM  could  enable  over  terahertz  communication  links.  In 
addition,  the  ability  of  SSH  to  couple  any  incident  beam  direction  to  any  outgoing  beam 
direction  depending  on  the  incident  beams  frequency  may  be  important  in  SSH  N  x  N  routing 
applications.  The  concept  is  described  in  the  proceedings  of  DEPOM’96  and  CLEO’97.  A 
manuscript  describing  a  demonstration  of  the  concept  is  being  submitted  to  Optics  Letters. 
Coherent  Transient  Simulator 

Several  of  the  work  done  under  this  program  have  relied  on  a  computer  program  which 
simulates  the  optical  coherent  transient  signal  produced  from  an  arbitrary  sequence  of  input 
excitation  pulses.  As  the  code  for  the  simulator  has  not  been  published  elsewhere,  it  is 
contained  in  appendix  A. 

The  versatile  simulator,  begun  under  AFOSR  contract  F49620-91-C-0088  and  improved 
under  that  current  grant,  integrates  the  optical  Bloch  equations  and  follows  the  evolution  of  the 
density  matrix  elements  during  and  between  excitation  pulses.  After  the  final  excitation  pulse, 
the  density  matrix  elements  contain  the  information  necessary  to  calculate  the  resultant  coherent 
transient  response  of  the  absorbing  medium.  Since  the  Bloch  equations  take  into  account 
nonlinear  effects,  the  program  is  useful  in  studying  the  effect  nonlinear  excitation  pulses.  The 
program  currently  runs  a  Gateway  486-33MHz  machine,  but  is  written  in  C  using  the  standard 
I/O  interface  and  is  thus  directly  portable  to  workstations  or  mainframe  computers  and  could 
be  incorporated  in  the  core  of  channel  modeling  program  used  to  evaluate  the  performance  of 
memory  and  processing  systems. 

The  coherent  transient  simulation  program  has  been  optimized  for  speed,  flexibility,  and 
input  and  output  capabilities.  The  simulator  accepts  as  input  an  arbitrary  number  of  excitation 
pulses,  including  chirped  pulses  and  arbitrarily  long  phase-,  amplitude-,  or  frequency-encoded 
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data  pulses.  The  exact  analytical  transformation  of  the  density  matrix  elements  resulting  from 
excitation  by  a  square  pulse  of  arbitrary  duration,  amplitude,  center  frequency,  phase,  and 
timing  was  calculated  and  incorporated  into  the  program.  This  greatly  enhanced  the  programs 
speed  since  integration  of  the  matrix  elements  was  no  longer  required.  Chirped  pulses  of 
arbitrary  duration,  bandwidth,  amplitude,  center  frequency,  phase,  and  timing  are  calculated 
by  integration  of  the  optical  Bloch  equations.  The  response  of  the  medium  to  other  types  of 
excitation  pulses  of  arbitrary  temporal  structure  can  be  calculated  in  the  same  way. 

Input  capabilities  and  calculation  flexibility  include  options  for  32-bit  data  pulses,  gating 
pulses,  sudden  losses  in  coherence,  and  homogeneous  decay.  Data  pulses  up  to  32-bits  long 
can  be  entered  as  a  single  integer  and  can  be  amplitude  or  phase  encoded,  return  to  zero  (RZ) 
or  non-return  to  zero  (NRZ).  Data  pulses  longer  than  32-bits  are  achieved  by  repeating  this 
option  indefinitely.  The  gating  pulse  option  simulates  the  effects  of  a  gating  pulse  by 
eliminating  those  atoms  in  the  upper  state  at  the  time  of  the  pulse,  thus  permanently  modifying 
the  inhomogeneous  profile.  The  resultant  modulations  in  the  inhomogeneous  profiles  will  alter 
the  effect  of  excitation  pulse  nonlinearities  compared  to  an  unmodulated  inhomogeneous 
profile.  The  ability  to  create  a  sudden  loss  of  coherence  aids  in  simulating  absorbers  with 
homogeneous  decay  times  much  shorter  than  the  upper  state  lifetimes  and  in  simulating  the 
rejection  of  spurious  outputs  via  spatial  isolation  (since  the  program  does  not  do  spatial 
integrations).  Homogeneous  decay  can  be  optionally  included  in  the  standard  Bloch  equations 
for  chirped  or  square  pulses.  The  output  signals  resulting  that  would  be  observed  with  direct 
detection  and  coherent  detection  are  calculated,  as  well  as  the  resultant  holebuming  spectra. 

The  graphical  interface  allows  quick  evaluation  of  the  results.  The  graphical  interface  is 
achieved  via  a  different  program  (Microsoft  Excel)  that  is  simultaneously  resident  with  the 
coherent  transient  simulator.  This  separation  maintains  the  quick  portability  of  the  simulator  to 
other  machines. 
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APPENDIXA 


CODE  FOR  COHERENT  TRANSIENT  SIMULATOR 

Some  of  the  following  routines  were  adapted  from  those  published  in  the  book  Numerical  Recipes  in  C  by  W 
H.  Press,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W.  T.  Vetterling,  Cambridge  Unviersity  Press,  1988.  Not 
included  here  are  routines  that  are  extensively  identical  to  those  found  in  Numerical  Recipes  in  C  and  one 
should  refer  to  this  book  for  the  code  for  these  routines. 

BLOCH.C 

/***  This  program  integrates  the  optical  Bloch  equations  for  a  two-level 
system  where  2(PI)V(t)/h  =  omegaR(t-tstart)*cos(omegaO*t  + 

(chirp/2)*(t-tstart)"'2  +  phiO  +  PM(t-tstart) ) 

=  omegaR(t-tstart)*cos(omegalO*t  +  phi(t)), 

whCTe 

detune  =  omega  10-omegac 

detuneL=  omegaO-omegac 

detuneAL  =  omega  10  -  omegaO  =  detune  -  detuneL 

phi(t)  =  A  +  B*(t-tstart)  +  C*(t-tstart)''2  +  PM(t-tstart) 

A  =  phiO  -  detuneAL*tstart 

B  =  -detuneAL 

C  =  chirp/2 

omega(t)=  omegaO  +  chirp*(t-tstart)  +  d/dt(PM(t-tstart)) 

=  instantaneous  frequency. 

Density  matrix  elements: 

rho00=Bloch[l],  rho01real=Bloch[2],  rho01imag=Bloch[3],population  density=Bloch[0] 

Pulse  parameters  that  can  be  changed  during  runtime  in  set_values()  are: 
omegaR,  chirp,  phiO,detuneL,tstart,tstop. 

detuneL  is  the  laser  detuning  from  center  of  the  inhomogeneous  profile, 
detune  is  the  atoms  detuning  from  the  center  of  Inh.  profile  and  varies 
from  detune  1  to  detuneN  where  Ndetune  is  the  number  of  detunings. 

Two  functions  are  needed  in  derivs: 
omegaR(t)=>omegaRt 
phi(t)=>phit. 

The  pro^am  integrate  the  QBE  for  Npuise  pulses  for  Ndetune  different 
detunings.  The  pulse  parameters  can  be  independently  set.  The  data 
is  appended  to  the  specified  output  file  after  each  pulse.  FFT  information  is 
stored  in  file  graphit.dat  for  graphing  with  IGOR. 

Project  includes: 

bloch.c,  complex.c,  OBEJnt.c,OBE_bs.c,OBE_rk.c,OBE_sq.c,OBE_lft.c, 
gasdev.c,  ranl.c,ran2.c,  and  nrutil.c 

#define  maxNdetune  4096+1  /***  maximum  number  of  detunings  ***/ 

#define  maxNN2  (2*(maxNdetune-l)) 

#define  maxNpulse  10  /**♦  maximum  number  of  pulses  ***/ 

#define  TINY  1.0e~30  /***  definition  of  a  tiny  value  ***/ 

#defineSQR(a)((a)*(a)) 

#define  PI  3.14159265358979 
#defme  PI2  6.28318530717958 
#include  <stdHb.h> 

#include  <conio.h> 

#include  <stdio.h> 

#include  <io.h> 

^include  <math.h> 

^include  "complex.h" 

^include  "bloch.h" 

#include  "nrutil.h” 

#include  "stdio  x.h” 
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FELE  *fp_save,*fp_read,*^_HB; 

float  huge  Bloch_matrix[l+maxNdetune][NVAR+l], 

float  CTsignal[maxNN2+l];  /***  defining  declaration  ***/ 

/***  The  following  values  can  be  modified  in  set_vaiues  routine  ***/ 
float  BIochO[NVAR+l]; 

double  tstart[maxNpulse-»- 1  ],tstop[maxNpulse+ 1  ], chirp BW[maxNpulse+ 1  ] ; 

double  omegaR[maxNpulse+l  ],phiO[maxNpuIse+l  ],detuneL[maxNpulse+ 1  ]; 

double  detune  1  ,detuneN,ddetune,phiD,eps,area,chirp,gamnia; 

double  noise,inhomo_HW; 

int  Ndetune,arun,NpulseJp,NN,NNpower; 

int  calc_mode[maxNpulse+l],data_flag[maxNpulse+l]; 

long  unsigned  int  data_vaiue[maxNpulse+l]; 

charsave_filename[20],read_filename[20]="",values_filename[20]-"”, 

charCT_^filename[20]=‘'graphit.t",HB_filename[20]=*'graphit.f’; 

int  save_flag,read_number, printout  flag; 

void  main(void) 

{  . 

int  brk,idetune,bit_valuej_sub; 
long  unsigned  int  mask; 
double  tau; 

/*  console_options.pause_atexit^O; 

console_options.nrows=37;  */ 

/**  Quick  C  for  Windows  additions:  ^*1 

/*_wabout('’Opticai  Bloch  Equation  SoIver\n\tby\nW.  R.  Babbitt*');  */ 

/*  *  end  of  Quick  C  Additions  *  *  / 

/*  BIoch_matrix=matrix(  1  ,maxNdetune,0,NVAR); 

CTsignal=vector(l,maxNN2);  */ 

/***  set  initial  v^ues  ***/ 
init_values(); 

/***  set  parameters  for  first  run  ***/ 
brk  =  set_values(); 
while((brk  =  0)) 

{ 

/***  Open  save  file  ***/ 
if{save_flag  =  I) 

^_save  =  fopen(save_filename,"a")  )=  NULL) 

printf("save_file  %s  can  not  be  opened\n",save  filename)' 
brk-1; 

save^flag  =  0; 

} 

/***  Initialize  Bloch__matrix  ***/ 
brk  +=  init_bIoch(); 

/***  print  run  parameters  to  screen  and  file 
prin4”%6s  %6s  %6s  %6s  %8s  %4s  %7s  %7s  %7s\n", 

''B10[l]","B10[2]","BI0[3]"/'phiD"/*eps”,"arun",*'Ndetune","detunei",'’detuneN")' 
printf("%6.3f %6.3f  %6.3f  %6.3f  %8g  %4d  %7d  %7.3f  %7.3f\n", 
BlochO[l],BlochO[2],BlochO[3],(phiD/PI),eps,arun,Ndetune,(detunel/PI2) 
(detuneNyTI2)); 

printf(”%6s=  %d\n’V'Npulse",Npulse); 
if  ((save_flag  =  l)&&(brk=0)) 

( 

fprintf(fp__save,"\n  \n"); 

f^rintf(f^_save,  "%6s\t%6s\t%6s\t%6s\t%8s\t%4s\t%7s\t%7s\t%7s\t%7s\n",  "BlOr  11" 
"B10[2]'V’B10[3]’V'phiD","eps","arun","Ndetune","detuner',"detuneN" 

'Tnho_HW"); 

fprintf(fp_save;'%6Jfit%63flt%6Jfit%6Jflt%8g\t%4d\t%7d\t%7Af\t%7.4f\t%7.4An", 

BlochO[l],BlochO[2],BIochO[3],(phiD/PI),eps,arun,Ndetune,(detunel/PI2), 
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(detuneN/PI2),(inhomo_HW/PI2)); 

fprmtfl;fp_save,"%6s\n","Npu]se"); 

4}rintf(4)_save,"%6d\n",Npulse); 

/***  integrate  pulses  ***! 

for(jp=l;((jp  <=  Npulse)&&(brk=0));  jp++) 

if  (data_flag|jp]==0)  /*♦*  Single  pulse  ♦**/ 

brk=integrate_pulse(Bloch_matrix,fp_save,jp,calc_mode[jp],eps,omegaR|jp] 
tstart[jp],tstop[jp],detuneL[jp],phiOOp],chirpBW[jp],gamnia, 
detune  1  ,ddetune,Ndetune,printout _flag,save_flag); 


else  /***  Data  stream  (multiple  subpulses)  ***/ 
for(j_sub=0,mask=(long  unsigned)  1  ;((mask<=data_value[jp])&&(brk=0));j_sub-H-,mask«=  1 ) 


tau=tstopOp]-tstart(jp]; 
bit_value=(mask&data_value[jp])?  1:0; 

printf("Pulse#%d  is  a  Data  Stream  data_value=%lu  mask=%lu  bit_value=%ld\n" 
jp,data_value[jp],mask,bit_value); 
switch(data_flag|jp]) 

{ 


case  1 ; 


brk=integrate_pulse(Bloch_matrix,fp_save,jp,2,eps, 
(omegaR|jp]*(double)bit_value),{tstart[jp]+(double)j_sub*tau), 
(tstart|jp]+(double)(j_sub+l)*tau),detuneL[jp],phiO[jp],cliirpBW[jp], 
0 . 0,  detune  1 ,  ddetune,  Ndetune,  0, 0); 
break; 


case  2: 

{ 

brk=mtegrate_pulse(Bloch_matrix,fp_save,jp,2,eps,omegaR[jp], 
(tstart[jp]+(double)j_sub*tau),(tstartOp]+(double)(j_sub+l)*tau), 
detuneL[jp],(phiO[jp]+PI*(double)bit_value).chirpBW[jp],0.0, 
detune  1  ,ddetune,Ndetune,0,0); 
break; 

} 

case  3: 

{ 

brk=integrate_pulse{Bloch_matrix,^_saveJp,2,eps, 
(omegaR[jp]*(double)bit_value),(tstart[jp]+(double)j_sub*tau), 
(tstartljp]+(double)(j_sub+ 1  )*tau),detuneL|jp], 
(phiOOp]+0.5*PI*(double)j_sub),chirpBW|jp],0.0, 
detunel,ddetune,Ndetune,0,0); 
break; 


case  4. 

{ 

brk=integratejDulse(Bloch_matrix,fp_save,jp,2,eps,omegaR[jp]* 
exp(((tstartOp]+(double)j_sub*2.0*tau)-20)*gamma*2.0), 
(tstart[jp]-^(double)j_sub*2.0*tau), 
(tstartOp]+((double)j_sub*2.0+1.0)*tau), 
detuneLOp],(phiO[jp]+PI*(double)bit_value),chirpBW[jp],0.0, 
detune  1  ,ddetune,Ndetune,  0,0); 
for(idetune=l;idetune  <=  Ndetune;  idetune-H-){ 

Bloch_matrix[idetune][2]*=exp(-2.0*tau*gamma); 
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B loch_niatrix[idetu ne]  [3  ]  *  =exp( -2 . 0* t au *  gamma) ; 

i 

break; 

} 

case  5: 

{ 

brk^integrate_pulse(BIoch_matrix,fp_save,jp,2,eps, 
(omegaR[jp]*(double)bit_value),(tstart[jp]+(double)j_sub*2.0*tau), 
(tstartljp]+((double)j_sub^2.0-M.0)*tau),detuneL[jp],phi0[jp], 
chirpBW[jp],0.0,detunel,ddetune,Ndetune,0,0), 
break; 
case  6: 

{ 

brk=integrate_pulse(Bloch_matrix,fp_save,jp,2,eps,omegaR[}p]* 

exp(-((tstart[jp]+(double)j_sub*2.0*tau)“35)*gamma*2.0), 

(tstartOp]+(double)j_sub*2.0*tau), 

(tstart  [jp]+((double)j_sub  *2 , 0+ 1 . 0 )  *  t  au), 
detuneL[jp],(phiO|jp]+PI*(doubie)bit__value),chirpBW[Jp],0.0, 
detune  1  ,ddetune,Ndetune,0, 0), 
for(idetune=  1  ;idetune  <=  Ndetune;  idetune^)  { 

Bloch_matrix[idetune][2]*=exp(-2.0*tau*gamma); 

Bloch_matrix[idetune][3]*=exp(-2.0*tau*gamma); 

} 

break; 

} 

case  7: 

{ 

brk=inte^ate_pulse(Bloch_matrix,fp_save,jp,2,eps,omegaR[jp], 
(tstart  [jp]+(double)j_sub*2.0*tau), 

(tstart  [jp]+((double)j_sub*2.0+L0)*tau), 
detuneL[jp],(phiO[jp]+PI*(double)bit_value),chirpBW[jp],0.0, 
detune  l,ddetune,Ndetune,0,0); 
for(idetune=l;idetune  <=  Ndetune;  idetune++){ 

Bloch_matrix[idetune][2]*=exp(-2.0*tau*gamma); 

Bloch_matrix[idetune][3]*=exp(-2.0*tau*gamma) 

} 

break; 

} 

case  8: 

{ 

brk=integratejulse(Bloch_matrix,fp_saveJp,2,eps,omegaR[jp]* 
exp(((tstartljp]+(doubIe)j_sub*2.0*tau)-20)*gamma*2.0X 
(tstart  [jp]+(double)j_sub*2.0*tau), 

(tstart  [jp]+((double)j_sub  *2 . 0+1.0)  *tau), 
detuneL^p],(phiO[jp]+PI*(double)bit_value),chirpBW[jp],0.0, 
detunel,ddetune,Ndetune,0,0); 
for(idetune=l;idetune  <=  Ndetune;  idetune++){ 

BIoch_matrix[idetune][2]*=exp(-2.0*tau*gamma); 
Bloch_matrix[idetune][3  ]  *=exp(-2. 0*tau*  gamma); 

} 

break; 

} 


}  /***  end  of  switch  ***/ 

}  /***  end  of  for(j_sub)  ***/ 

/*printf("\n");*/ 
if  (printout_flag  =  1) 

prmtf("%4s  %7s  %9s  %10s  %10s  %9s\n","idet'\ 
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"detune'',"Bloch[l]","Bloch[2]","Bloch[3]","Bloch[0]"); 

if  (saveflag  ==  1) 

{ 

area=  omegaR[jp]*tau; 

chirp=(fabs(chirpB  W[jp])>TINY)?(chirpBW[jp]/tau):0 . 0; 
fprinlf(fp_save,"%6s\t%4s\t%4s\t%10s\n","pulse#",”calc","data","data_value"); 
fi)rintf(fj)_save,"%6d\t%4d\t%4d\t%  1 01u\n" jp,2,data_flag[jp],data_value0p]);  ’ 
fprintf(fi)_save, "%  1 2s\t%8s\t%8s\t%  1 0s\t%7s\t%5s\t%7s\t%  1  Os  \n", 
"omegaR”,"tstart","tstop","area","detuneL","phiO","chirpBW","chirp"); 
f)rintf(fp_save,"%12.6g\t%8.3f\t%8.3f\t%10.3g\t%7.3f\t%5.3f\t%7.3ftt%io.3g\n”, 
(omegaR0p]/PI2),tstart[jp],tstop[jp],(area/PI),(detuneL[jp]/PI2), 
(phiO|jp]/PI),(chirpBW[jp]/PI2),(chirp/PI2)); 

^rintf(fp_save,  "%4s\t%7s\t%9s\t%  1  Os\t%  1 0s\t%9s\t%4s\t%4s\n", 
"idet","detune'',"Bloch[l]","Bioch[2]’',’'Bloch[3]”,"Bloch[0]",''junk",”junk"); 

for(idetune=l;((idetune  <=  Ndetune)&&(brk==0));  idetune++) 
if  (printoutfiag  =  1) 

printf("%4d%7,4f%9.7f%]0.7f%10.7f%9.7f\n'',idetune, 

((detune  1  +(double)(idetune-  l)*ddetune)/PI2),Bloch_matrix[idetune][  1  ], 
Bloch_matrix[idetune][2],BIoch_matrix[idetune][3]); 
if  ((save  flag  ==  1  )&(printout  flag  !=  2)) 

fi)rintf(fp_save,''%d\t%.4g\t%g\t%g\t%g\t%g\t%d\t%d\n", 
idetune,((detune  1  +( double)(idetune- 1  )*ddetune)/PI2), 

Bloch_matrix[idetune][  1  ],Bloch_matrix[idetune][2], 
Bloch_matrix[idetune][3],Bloch_niatrix[idetune][0],0,0); 

}/***endofifelse***/ 

}  !***  end  of  for(jp)  ***! 

I***  Store  Holebuming  Spectrum  out  for  graphit  to  graph  ***/ 
ifl:brk=0) 

{  , 

if((  fp_HB  =  fopen(HB_filename,"w")  )=  NULL) 

printf("HB_File  %s  can  not  be  opened", HB  filename)' 
brk=l; 

} 

else 

{ 

for(idetune=l;(idetune  <=  Ndetune);  idetune++) 

fprintf(fp_HB,  "%g\t%g\n", 

((detunel+(double)(idetune-l)*ddetune)/PI2), 

Bloch_matrix[idetune][l]*Bloch_matrix[idetune][0]), 

fclose(fp_HB); 

} 

} 

/*♦*  Do  FFT  if  no  break  and  Ndetune  >  16  and  frequency  center  is  0  GHz  ***/ 
if  ((brk  =  0)&&(Ndetune>16)&&(SQR(detuneN+detunel)<TINY)) 
brk=FFT_output(Bloch_matrix,CTsignal,detunel,detuneN,inhomo_HW,noise, 

phiD,NN,CT_filename,0);  /*gamma)  Moe  12/19/96;*/ 
if  (save  flag  =  1 ) 
fclose(^_save); 
if  ((arun  <=  l)||(brk!=0)) 

{  . 

printf(  "Done*  **’*'*************************\jQ‘»y 
/*  printf("Done%c%c%c\n";\7^^7^^7);  */ 
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/ 


brk=set_values();  /♦**  set  parameters  for  next  run  ***/ 

} 

else 

brk=autonan()y***  perform  arun  commands  ***/ 

/***  Close  output  file  after  each  run  ***/ 

}  /***  end  while(brk  =  0)  ***/ 

/*  free_matrix(B!och_matrix,  1  ,maxNdetune,0,NVAR); 
fi-ee_vector(CTsignal,l,maxNN2);  */ 

int  autorun(void) 

{ 

int  brk=0; 

jp=i; 

arun-;  /***  decrement  number  of  automatic  runs  ***/ 

/*omegaR[l]  *=0.1; 
area=  omegaR[  1  ]  *(tstop[  I  ]-tstart[  1  ]);*/ 
retum(brk); 

} 

int  init_bloch(void) 

{ 

int  brk=0,idetune,ivar,iline,  cbuff,  i  occur; 
float  dummy  1 ,  dummy2,  dummy3 ,  dummyO 
charibufiI256]; 

/***  Initialize  Bloch_matrix  ***/ 
fbr(idetune=  1  ;idetune  <=  Ndetune;  idetune++) 
for(ivar=0;ivar  <=  NVAR;  ivar-H-) 

Bloch_matrix[idetune][ivar]=BlochO[ivar]; 

if(read_number>0) 

{ 

if  (save  flag  =  1) 
fclose(fp_save); 

ifl:(fp_read  =  fopen(read_filename,"r")  )!= NULL) 

printf("Searching  file  %s  for  %d  occurrence  of  idetune  string.\n", 
read_filename,read_number); 
for(i_occur=  1  ;i_occur<=read_number,i_occur-H-) 
for(iline=  1  ;((fgets(ibufi;  1 00,^_read)!=NULL)&&!(sscanf(ibufF,"idet%c", 

&cbuff)));iline-H-); 

if(!(sscanf(ibufl;"idet%c",&cbuff))) 

{ 

if(i_occur=l) 

printf("The  string  idetune  was  not  found  in  file  %s  \n",read_filename); 
else 

printlC'Only  %d  occurrences  of  idetune  string  were  found\n",i_occur); 
read_number=0; 
brk=l; 

} 

else 

{  . 

printf( "Reading  file  %s  starting  at  line  %d\n",read  filename  iline+1)' 
fgets(ib^,255,fp_read); 

sscanflibufiF,  "%d  %lf  %f  %f  %f  %f ',&idetune,&detune  1  ,&dummy  1  ,&dummy2, 
&dummy3,&dummy0); 

Bloch_matrix[  1][  1  ]=dummy  1 ; 

Bloch_matrix[l][2]=dummy2; 

Bloch_matrix[  1]  [3]=dummy3 ; 

Bloch_matrix[  1  ]  [0]=dummy0 , 
detunel  *=  PI2; 
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while  ( (fgets(ibuffJOO,fp_^read)!=NULL)  &&(sscanf(ibuff,”%d”,&idetune)) ) 
sscanf(ibulf;’%d%If%f%f%f%f\&idetune,&detuneN,&dummyl,&dummy2, 

&dummy3  ,&dummyO); 
if  (idetune!=- 1 ) 

f 

I 

Blochmatrix  [idetune]  [  1  ]=dummy  1 ; 

Blochmatrix  [idetune]  [2]=dummy2, 

Bloch_matrix[idetune][3]=dummy3; 

Bloch_matrix[idetune][0]=dummyO' 

} 

} 

if  (idetune=- 1 ) 

{ 

printf("Error  in  reading  read  file,  idetune  =  -1.  Cancelling  read_file.\n"); 

detuneN=ddetune*((doubleXNdetune- 1  )>+detune  1 ; 

readnumber^O; 

brk=l; 

} 

else 

{ 

detuneN  *=  PI2; 

Ndetune=idetune; 

NN=Ndetune-l; 

if(Ndetune  !=  l)ddetune=(detuneN-detunel)/((double)(Ndetune-l)); 
printf("Ndetune=%7d  detunel/PI2=%7.3f  detuneN/PI2=%7.3f\n",Ndetune 
(detunel/PI2),(detuneN/PI2)); 
if  (printout  flag  =  1 ) 

{ 

printfl:"%4s  %7s  %9s  %10s  %10s  %9s\n","idet", 
"detune","Bloch[l]”,"Bloch[2]","Bloch[3]","Bloch[0]"); 
for(idetune=l;idetune  <=  Ndetune;  idetune-H-) 
printf("%4d%7.4f%9.7f%10.7f%10.7f%9.7f\n",idetune, 

((detune  1  +(double)(idetune- 1  )*ddetune)/PI2),Bloch_niatrix[idetune][  1  ], 
Bloch_matrix[idetune][2],Bloch_matrix[idetune][3], 
Bloch_matrix[idetune]  [0]); 

} 

} 

} 

fclose(fpread); 

} 

else 

{ 

printf("Read_File  %s  can  not  be  opened ",read  filename); 

read_number=0; 

brk=l; 

} 

if  (save_flag  =  1) 

f 

V 

ifl;(  fpsave  =  fopen(save_filename,"a")  )=  NULL) 

printf("save_file  %s  can  not  be  reopened",save  filename)' 
brk=l; 

save_flag  =  0; 

} 

else 

{  . 

rf(read_number>0) 
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{ 

fprintf(fp_save,"Read_fiIe=  %s  Occurrence^  %d  Starting  line=  %d\n", 
read_filename,read_number,iline+ 1 ); 
fprintf(fp_save,  "%4s\t%7s\t%9s\t%  1  Os\t%  1  Os\t%  1  Os\n", 
”idet^"detune^"Bloch[l]^”Bloch[2]^"Bloch[3]^”Bloch[0]”); 
for(idetune=l;idetune  <=  Ndetune;  idetune-H-) 
fprintfl^fp^save,  ”%4d\t%7.4f\t%9. 7f\t%  1 07f\t%  1 0. 7f\t%  1 0.7f\t%d\t%d\n", 
idetune,((detune  1  +(doubIe)(idetune- 1  )*ddetune)/PI2), 
Bloch_matrix[idetune][  1  ],Bloch_matrix[idetune][2], 
Bloch_matrix[idetune][3]31och_niatrix[idetune][olo,0); 

} 

} 

} 

} 

retum(brk); 

}  . 

void  init_values(void) 

{ 

anin=0;  /***  arun  =  number  of  times  program  is  to  run  automatically  ***/ 

printout_flag=0; 

read_number=0; 

read_filename[0]=’\0’; 

phiD=0; 

NNpower=l; 

NN  =  1  «  NNpower; 

Ndetune=NN+l;  /*♦*  number  of  detunings  to  evaluate  =  1+  2^NNpower  ***/ 

detunel  =  -2.000*PI2;  /***  first  detuning  in  Grad/sec  ***/ 

detuneN  =  2.000*PI2;  /***  last  detuning  in  Grad/sec  ***/ 

inhomo_HW  =  0.0; 

noise  =  0.0; 

if(Ndetune  !=  1)  ddetune=(detuneN“detunel)/((doubleXNdetune-l)); 

Bloch0[l]=1.0;  /***  starting  value  of  rhoOO  ***/ 

BIoch0[2]=0.0;  /***  starting  value  of  real  part  of  rhoOl  ***/ 

Bloch0[3]=0.0;  /***  starting  value  of  imag  part  of  rhoOl  ***/ 

Bloch0[0]=1.0;  /***  starting  value  of  population  ***/ 
eps=1.0e-9;  /***  Set  error  parrameters  ***/ 

Npulse=l ;  /*♦*  number  of  pulses  to  integrate  ***/ 

for(jp=l;jp  <=  maxNpulse;  jp-H-) 

{ 

data  flag[jp]=0; 

data_value[jp]=(long  unsigned)0; 

tstart|jp]=5.0;  /***  tstart,tstop  in  nanoseconds  ***/ 

tstop[jp]=15.0; 

omegaR[jp]=0.025*PI2;  /***  omegaR  in  gigarads/sec  ***/ 

phi0[jp]=  0.0*PI;  /***  phiO  in  radians  ***/ 

chirpBW[jp]=  0.000*PI2;  /***  chirpBW  in  gjgarads/sec  ***/ 

detuneL[jp]=  0.000*PI2;  /***  detune  in  gigarads/sec  ***/ 

calc_mode[jp]=2;  /**♦  o  for  Runge-Kutta,  1  for  Bulirsch-Stoer,  2  for  analytic  ***/ 

tstart[2]=45.0; 

tstop[2]=55.0; 

omegaR[2]=0.05*PI2; 

gamma=0.0; 

jp=l ;  /***  current  pulse  number  ***/ 

chirp=chirpBWOp]/(tstopOp]  -  tstartQp]);  /*♦*  chirp  in  gigarads/sec/ns  ***/ 
area=  omegaR[jp]*(tstop[jp]-tstart[jp]);  /***  pulse  area  of  current  jjp]  pulse  ***/ 

int  set_vaiues(void) 
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{ ^ 

int  brk=0,jjpj_sub,jp_insert,bit_vaIue; 

unsigned  long  int  mask; 

char  choice,ibufff256]; 

int  nchoice; 

char  choicel,choice2; 

char  schoice[80]; 

double  ytemp; 

arun=0; 
choice=  '\0’; 
while  (choice  !=  V) 

{ 

area=omegaRljp]  *(tstop  [jp]-tstart  Op]); 

chirp=(fabs(chirpBW[jp]>TINY))?(chirpBWOp]/(tstop[jp]  -  tstart[jp])):0.0; 
if((choice  =  V')|l(choice='V)) 

{ 

printf("%6s  %6s  %6s  %6s  %8s  %4s  %7s  %7s  %7s  %7s\n","B0[l]","B0[2]", 
"B0[3]'V'phiD","eps","arun","Ndetune","detunel","detuneN'',"Inho_HW"); 

printf("%6s  %6s  %6s  %6s  %8s  %4s  %7s  %7s  %7s  %7s\n" , 
"rho00","rho01r","rho01i","{/PI}"," ",”{GHz}","{’gHz}","{GHz}"); 
printf("%6.3f  %6,3f  %6.3f  %6.3f  %8.  le  %4d  %7d  %7,3f '%7.3f  %7.3f\n", 
BlochO[l],BlochO[2],BlochO[3],(phiD/PI),eps,arun,Ndetune,(detunel/PI2) 
(detuneN/PI2),(inhomo_HW/PI2)); 
printf("%6s  %6s  %6s  %6s  %8s  %4s  %7s  %7s  %7s  %7s\n", 

"Z,z","X,x","Y,y","P","E,e","U,u'',"N,n","D,d'',"D,d","Ww")- 

printf("  \n"); 

printf(’’%4s  %4s  %4s  %4s  %10s  %4s  %4s%16s%16s\n","Npls","pls#","mode","data", 
"data_value","pmt","rd#","read_filenanie","save_filename"); 
prmtfl:"%4s  %4s  %4s  %4s  %  I  Os  %4s  %4s%  1 6s®/ol  6s\n", 

printf("®/;4d  %4d  %4d  %4d  ®/ll6lu’%4d  %4d%16s%16s\n",Npulse,jp,calc_mode[jp], 
data JlagUp],  data_valueljp],printout_flag,read_number,read  filename, 
save_filename); 

printf("%4s  %4s  %4s  %4s  %10s  %4s  %4s%16s%16s\n",”K  k"  "J  i"  "Mm"  "HOhl"  "H  h" 
"F,f',"r","R->param,r->rho","S,s")  .  ,  ,  , 

printf("\n"); 

printf("%8s  %8s  %8s  %IOs  %7s  %5s  %7s  %10s  %8s\n", 

"omegaR","tstart","tstop","area","detuneL","phiO","chirpBW","chirp","Homo  LW”) 

printf("%8s  %8s  %8s  %  1  Os  %7s  %5s  %7s  %  1  Os  %8s\n", 

"{GHz}","{nsec}","{nsec}","{/PI}","{GHz}","{/PI}","{GHz}","{GHz/ns}"  "(GHz)") 
prmtf("%8.6g  %8.3f  %8.3f  %10.3g  %7.3f  %5.3f  %7.3f  %10.3g  %8e\n", 
(omegaR[jp]/PI2),tstart|jp],tstop[jp],(area/PI),(detuneL[jp]/PI2), 
(phiO[jp]/PI),(chirpBW0p]/PI2),(chirp/PI2),(gamma/PI)); 
printf("%8s  %8s  %8s  %10s  %7s  %5s  %7s  %10s  %8s\n", 

"0,o","T,t","T,t","A->t,a->o","L,l","p","B,b->L","C->t,c->b"  "G  s")- 

if(choice='V) 

{ 

printf("\nRecap  of  all  %d  pulses\n",Npulse); 
for(jjp=  1  Jjp<=Npulse;jjp++) 

{ 

printf{"%4d  %4d  %4d  %4d  %101u\n",Npulse,ijp,calc_mode[ijp],data_flag[jyp], 
datavalueQjp]); 

printfl:"%12.6g%8.3f%8.3f%10.3g  %7.3f%5.3f %7.3f %10.3g\n", 
(omegaR[ijp]/PI2),tstart[ijp],tstop[ijp], 
((omegaR[ijp]*(tstop[ijp]-tstart[ijp]))/PI),(detuneL[ijp]/PI2X 

(phiO[ijp]M),(chirpBW0jp]/PI2), 

(((fabs(chirpBW[ijp])>TINY)?(chirpBW{ijp]/(tstop[jjp]  -tstart[ijp]));0.0)/PI2)); 
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} 

} 

printfi:"  \n"); 

}  . 

printfi;"Enter  choice  (v=variable  Iist,q=quit,i=integrate)(currently,  jp=%d):  "jp); 

choice  =  (char)getche(); 

printfl;"\n"); 

switch(choice) 

{ 

case  'A': 
area  /=  PI; 
do 
{ 

printf("input  area  (in  units  of  PI)(>=0)  (modifies  tstop):  "); 
sscanf(gets(ibuli),"%lf',&area), 

}whi]e  (area<0.0); 
area  *=  PI; 

tstop[jp]=tstart[jp]  +  ((omegaR[jp]>TINY)?(area/omegaR[jp]):TINY), 
break; 
case  'a': 
area  /=  PI; 
do 

{  , 

printf("input  area  (in  units  of  PI)(>=0)  (modifies  rabi  frequency):  "); 
sscanl(gets(ibuff),"%lf',&area); 
jwhile  (area<0.0); 
area  *=  PI; 

omegaROp]  =  area/(tstop[jp]-tstart[jp]); 
break; 
case  'B'; 

chirpBWQp]  /=  PI2; 

printf("input  chirpBW  (in  gigahertz)(Doesn't  modify  detuneL): ") 

sscanf(gets(ibufi),"%lf',&chirpBWOp]); 

chirpBWQp]  ♦=  PI2; 

chirp=(fabs(chirpBW[jp]>TINY))?(chirpBW[jp]/(tstop[jp]  -  tstart[jp]));0.0; 
break; 
case  'b': 

chirpBWQp]  /=  PI2; 

printfC’input  chirpBW  (in  gigahertz)(also  modifies  detuneL)  ") 

sscanf(gets(ibufi),"%lf',&chirpBWOp]); 

chirpBW[jp]  *=  PI2; 

chirp=(fabs(chirpBWOp])>TINY)?(chirpBW[jp]/(tstop|jp]  -  tstartnp])):0.0; 

detuneL[jp]  =  ~0.5*chirpBW[jp]; 
break; 
case  'C’; 
chirp  /=  PI2; 
do 

printf("input  chirp  (in  gigahertz/nsec)(naodifies  tstopjjp]);  "); 
sscanf(gets(ibufi),  ”%If ’,&chirp); 
jwhile  (chirp=0.0); 
chirp  *=  PI2; 

tstopOp]  =  (chirpBW[jp]/chirp)+tstart[jp]; 
break; 
case  ‘c’: 
chirp  /=  PI2; 

printf("input  chirp  (in  gigahertz/nsec)(modifies  chirpBW[jp]):  "); 
sscanf(gets(ibuj5),  "%If ’,&chirp); 
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chirp  *=  PI2; 

chirpBWljp]=chirp*(tstop[jp]-tstart|jp]); 
break; 
case  'D'; 
case ’d': 

il(read_number=0) 

{ 

do 

{ 

detune!  /=  PI2; 
printf("input  detune  I:  "); 
sscanfl;gets(ibufi),  "%lf' ,&detune  1 ); 
detune!  *=PI2; 
detuneN  /=  PI2; 

printfC'input  detuneN  (must  be  >  detune!):  "); 
sscanf(gets(ibuff),  "%If ’,«&detuneN); 
detuneN  ♦=  PI2; 

}  while  (detune!  >=  detuneN); 
if  ((Ndetune>16)&&(SQR(detuneN+detunel  )<TINY)) 
printf("FFT  info:  deita_t  =  %6.3f  nsec  Maximum  t  =  %8.3f  nsec", 

PI2/(detuneN-detunel),((double)(NN-l)*PI2/(detuneN-detune!))); 

else 

printf("Can’t  set  detunings  when  read  number  =%d.\n",read  number) 
break; 
case  'E': 
case  'e': 

printfC’input  eps:  "); 
sscanf(gets(ibuflF),"%lf’,&eps); 
break; 
case  'F: 
case  'f: 
do 
{ 

printfCIFrequency  Information  Printout  Options:\nO)  File  Printout  Only\n"); 
printfC'l)  File  and  screen  printout\n2)  no  printout  to  file  or  screen\n"); 
printfC'input  desired  option:  "); 
sscanf(gets(ibuff),"%d",&printout_flag); 

}whiie  ((printout_flag!=0)&&(printout_flag!=l)&&(printout_flag!=2)); 

break; 
case  'G*: 
case  'g': 

if  (gamma=0.0) 

{ 

printf(” Currently  T2  is  infinite  and  gamma  =0\n’'); 
ytemp=1.0e9; 

} 

else 

{ 

ytemp  =  1.0/gamma; 

^  printf("Currently  T2  =  %lf  nsec  \n’', ytemp); 

printfC'input  homogeneous  decay  time,T2  (  >=Ie9  means  inf):  "); 
sscanf(gets(ibuff),  "%lf ',&ytemp); 
if  (ytemp>=I.0e9) 

{ 

gamma=0.0; 

printf(['T2  is  now  infinite  and  gamma  =0\n”); 
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else 


{ 

gamma  =  1.0/ytemp; 

printfT  T2  =  %lf  nsec,  gamma  =  %lf,  Homo  LW=  %lf  GHz\n",ytemp,gamma,(gamma/PI)); 

break; 
case  'H": 
data_flag(jp]=0; 

data_value[jp]=(long  unsigned)0; 
break; 
case  'h': 
do 
{ 

printfCmodulation  mode  options:  \n''); 

printf("0)Cancel  1)AM  NRZ  2)PM  3)AM  w/phase  rotation  4)PM  w/RZ  &  up  ramp  &  decay\n") 
printf("5)AM  RZ  6)  PM  w/RZ  &  down  ramp  &  decay  7)PM  w/RZ  w/decay:\n"); 
printf("input  modulation  mode:  "); 
sscanf(gets(ibuflE),"%d",&data_flag[jp]); 

}while((data _flagQp]<0)&&(data_flagOp]>7)); 
printf("input  data_value  (0  to  4294967295):  "); 
sscanf(gets(ibuff),"%lu",&data_vaiue[jp]); 

for(j_sub=0,mask=(longunsigned)l;((mask<=data_valueOp])&&(brk=0));j_sub++,mask«=l) 

bit_value=(mask&data_value(jp])?l:0; 
printfr"%  1  d",bit_value); 

} 

printf("  is  the  binary  sequence  for  %lu’',data_value[jp]); 
calc_modeljp]=2; 
break; 
case  T: 
do 
{ 

printf("Input  #  of  pulse  to  insert(>0)/deleted(<0)/escape(0):"); 
sscanf(gets(ibuff),  "%d",&jpjnsert); 
if  (abs(jp_insert)>Npulse) 
print^"abs(jp_insert)  >  Npulse  =  %d\n",Npulse); 
jwhile  (abs(jp_insert)>Npulse); 
if  Op_insert<0) 


for(ijp=-jp_insert;jjp  <  Npulse;  jjp-H-) 

{ 

data_flag[iijp]=data_flag[ijp+ 1  ]; 
data_value[j[jp]=data_value[ijp+ 1  ]; 
tstart[jjp]=tstart[ijp+ 1  ]; 
tstop[jjp]=tstopOjp+ 1  ]; 
omegaR[ijp]=omegaR[j[jp+ 1  ]; 
phi0[jjp]=phi01jjp+l]; 
chirpBWpjp]=chirpBW[jjp+ 1  ]; 
detuneL[ijp]=detuneL[ijp+l  ]; 
calc_modeQjp]=calc_modeQjp+l  ]; 

}  . 

printf("Pulse#  %d  to  pulse#  %d  have  been  shifted  down  one.\n", 
(-jpinsert),Npulse); 

printf("The  old  pulse#  %d  has  been  deleted.\n",(-jpjnsert)); 

Npulse—; 

break; 

} 
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if  (Npulse  maxNpulse) 

{ 

printf( "Npulse  =maxNpuIse=  %8d.  Can’t  insert  more  pulses.  \n", maxNpulse); 
break; 

) 

if  (jpjnsert>0) 


for(jjp=Npulse;jjp  >=  jp  insert;  ijp-) 

{ 

data JlagQjfrHl  ]=(iata Jlag[ijp]; 
data_value[jjp+ 1  ]=data_value[ijp]; 
tstart[jjp-Hl  ]=tstart[ijp]; 
tstop[jjp+l]=tstopOjp]; 
omegaR[jjp+ 1  ]-omegaR0jp]; 
phiO[jjp+l]=phiO[ijp]; 
chirpB  Wpjp+ 1  ]=chirpBW|jjp] ; 
detuneLQjp+1  ]=detuneL[j[jp]; 
calc_model]jp+ 1  ]=caic_mode[ijp]; 

}  . 

printf( "Pulse#  %d  to  pulse#  %d  have  been  shifted  up  one. \n"jp_insert, Npulse); 
printf(’’Pulse#  %d  is  a  duplicate  of  the  old  pulse#  %d.\n",jp_insert,jpjnsert); 
Npulse-H-; 
break; 

} 

break; 
case  'i': 

printf("Integrating  (press  b  to  break)"); 
break; 
case  7': 
case  'j': 

jp  =  Npulse  +1; 
while  (jp  >  Npulse) 

{ 

printf("input  pulse#  (=jp):  "); 

sscanf(gets(ibuff),  "%d",&jp); 

if  (jp  >  Npulse)  printf("Npulse=  %8d  \n",Npulse); 

} 

area=omegaR|jp]*(tstop[jp]-tstart[jp]); 

chirp=(fabs(chirpBWOp]>TINY))?{chirpBW[jp]/(tstop[jp]  -  tstartOp])):0.0; 
break; 
case  'K'; 
case  V: 

Npulse  =  maxNpulse  -^1; 

while  ((Npulse  >  maxNpuIse)||(Npulse<l)) 

{ 

printf("input  Npulse:  "); 
sscanft;gets(ibu^,  "%d",&NpuIse); 

if  (Npulse  >  maxNpulse)  printf("maxNpulse=  %8d  \n",maxNpuise); 
jp=i; 


break; 
case  'L’: 
case  T: 

detuneLjjp]  /=  PI2; 

printf("input  detuneL:  "); 

sscanf(gets(ibuff),"%lf',&detuneL[jp]); 

detuneL[jp]  *=  PI2; 

break; 
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case  'Nf ; 
case 'm': 

ijfi[data_flag|jp]=0) 

{ 

do 

{ 

printf("%s\n%s\n%s\n%s\n%s\n%s^\n%s\n%s\n’^ 

"0=Runge-Kutta’',"l=BuIirsch-Stoer  (constant  amp)", 

"2=Square  Pulse","3=Coherence  Loss","4=Gating  Pulse", 

”5=Bulirsch-Stoer  ^v^th  decay","6=Bulirsch-Stoer  (shaped  pulse)", 

"7=gating  and  noise"), 
printfC'input  caIc_mode:  "); 
sscanf(gets(ibuff),"%d",&calc__mode[jp]); 
jwhile  ((calc_modeOp]<0)|j(calc_mode[jp]>7)); 
if  (calc_mode[jp]=5) 

printf(" Currently  Homo  LW=  %lf  GHz  \n",gamma/PI); 
if((caIc_mode|jp]=3)||(calc_mode[jp]==4)||(calc_mode|jp]=7)) 

tstart[jp]=(jp>  1  )?tstop[jp- 1  ]  :0 . 0; 
tstop[jp]=tstart[jp]+0.00 1 ; 
omegaR[jp]=0.0; 
chirpBW[jp]=0.0; 

} 

} 

else 

printf("Can't  change  calc  mode  unless  data  flag=0  \n")' 
break; 
case  TsP: 
case  'n*: 

if(read_number=0) 

{ 

do 

{ 

if  ((NdetunO  1 6)&4&(SQR(detuneN+detune  1  )<TINY)) 
printf("NNpower=%d  Ndetune=%d  deitaj=%6.3fhsec  Maximum  t=-% 8.3 ftisec\n", 
NNpower,Ndetune,PI2/(detuneN-detune  1 ), 

((double)(NN- 1 )  *PI2/(detuneN-detune  1 ))); 
printf("input  NNpower  (Ndetune=2^NNpower  +1)"); 
sscanf(gets(ibuff),  ”%d",&NNpower); 

NN  =  1  «  NNpower; 

Ndetune=NN+l; 

if  ((Ndetune>maxNdetune)||(NNpower>  14)) 
printf("maxNdetune=  %8d  \n",maxNdetune); 
jwhile  ((Ndetune>maxNdetune)||(NNpower>14)); 
printfC’Ndetune  =  %d\n",Ndetune); 
if  ((Ndetune>  1 6)&&(SQR(detuneN+detune  1  )<TINY)) 
printff'TFT  info:  delta_t  =  %6.3f  nsec  Maximum  t  =  %8.3f  nsec\n", 

^  PI2/(detuneN-detune  1  ),((double)(NN- 1  )*PI2/(detuneN-detune  1 ))); 

else 

printfC'Can’t  set  Ndetune  when  read  number  =%d.\n",read  number)- 
break; 
case  'O': 
case  'o': 

omegaRUp]  !=  PI2; 
do 
{ 

printf("input  omegaR  (in  gigaghertz)(>=^0):  "); 
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sscanf(gets(ibuff),"%lf',&omegaR[jp]); 

}  while  (omegaR[jp]<0.0); 
omegaR[jp]  *=  PI2; 
area==omegaR|jp]*(tstop[jp]-tstart[jp]); 
break; 
case  ’P’: 
phiD  /=  PI; 

printf("input  phase  of  coherent  detection  phiD  (in  units  of  PI):  "); 
sscanf(gets(ibuff),"%lf',&phiD); 
phiD  *=  PI; 
break; 
case  ’p': 
phiO[jp]  !=  PI; 

printfC’input  phiO  (in  units  of  PI):  "); 
sscanI{gets(ibuf!)/'%lf',&phiO[jp]); 
phiO[jp]  *=  PI; 
break; 
case  'Q’: 
case  'q': 
choice-i'; 
brk=l; 
break; 
case  'R': 

brk=read_values(); 
choice- V; 
brk=0; 
break; 
case  'r‘: 

printfC’Input  file  from  which  to  read  in  Bloch  values  (CR=none):  "); 

gets(read_filename); 

if  (read_filename[0]  =  ’\0') 

{ 

printf(’'No  entry  (setting  read_number=0):\n  ”); 
read_number=0;  /***  No  output  file  ***/ 

else 

{ 

readnumber=0; 

printf("Input  occurance  number  of  desired  bloch  values  (read_number>0):  "); 
sscanf(gets(ibuff),  ’'%d",&read_number); 

} 

if(readnumber>0) 

printf("WiIl  read  file  =  %s  for  occurance  #%d  \n”,read_filename,read_number); 
else 

{ 

read_number=0; 

printf("Occurance  number  =0.  Cancelling  read_file.\n  "); 
read_filename[0]  =  ^0’; 

} 

break; 
case  'S’: 
case  's': 

printf("Input  the  save  file  name  (type  CR  for  none):\n”); 
gets(save_filename); 
if  (save_filename[0]  =  (char)NlJLL) 
save_f[ag=0;  /***  No  output  file  ***/ 
else 

{ 
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save_flag  =  1; 

printf(”output  file  =  %s  \n’',save_filename); 

} 

break; 
case  T'; 
case ’t’: 
do 

{ 

printf("input  tstart  (in  nanseconds):  "); 
sscanf(gets(ibuff),''%lf',&tstart[jp]); 
printfC’input  tstop:  "); 

sscanf(gets(ibuff),"%lf’,&tstopOp]); 

}  whiIe(tstart[jp]>=tstop[jp]); 
area=omegaR[jp]*(tstop|jp]-tstart[jp]); 

chirp=(fabs(chiipBW[jp])>TINY)?(chirpBWOp]/(tstopOp]  -  tstart[jp])):0.0; 

break; 
case  TT: 
case  V: 

printf("input  arun  (l=autorun):  "); 
sscanf(gets(ibuff),’'%d",&arun); 
break; 
case  W’: 

inhomo  HW  /=  PI2; 

printf("input  inhomo  HW  (0:flat  >0:Guassian  HWl/eM,<0:Lorentzian  HWHM)(GHz):  ”) 
sscanf(gets(ibuff),  ”%If  \&inhomo_HW); 
inhomo  HW  *=  PI2; 
if  (SQR(inhomo_HW)<TINY) 
printf("Flat  inhomogeneous  profile  \n"); 
if  (inhomo_HW>TIN  Y) 

printf("Guassian  inhomogeneous  profile  with  HWl/eM  =  %lf',inhomo  HW/PI2)* 
if  Onhomo_HW<(-TrNY)) 

printfC’Lorentzian  inhomogeneous  profile  with  HWHM  =  %lf',(-inhomo_HW/PI2)); 
break; 
case  'w': 
noise  *=  PI2; 

printf('^noise=%lf  per  GHz  .  Input  new  noise  per  GHz:  ", noise); 
sscanf(gets(ibuff),  ”%lf ',&noise); 
noise  /=  PI2; 
break; 
case  ‘X’: 
case  'x': 

printfl;"input  Bloch0[2]  (rhoOlr):  "); 
sscanf(gets(ibuff)/’%f',&BlochO[2]); 
break; 
case  ’Y’; 
case  'y': 

printfC’input  Bloch0[3]  (rhoOli):  "); 
sscanf(gets(ibuflO,”%f*,&BlochO[3]); 
break; 
case  ’Z’: 
case  'z': 

printf(’’input  BiochO[l]  (rhoOO):  "); 
sscanf(gets(ibufF),"%f ',&BlochO[  1  ]); 
break; 

} 

if(Ndetune  !=  l)ddetune=(detuneN-detunel)/((double)(Ndetune“l)); 
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^ “^)&&((calc_mode|jp]  !-5)&&(calc_mode[jp]  !=6)&&(calc_mode[jp]  != 1  )&&(caic  mode[jp] ! 

{ 

printf("calc_mode  must  equal  0  or  1  or  5  or  6  when  chirp  !=  0.  Setting  calc_mode=l\n''); 
calc_mode|jp]=l; 

) 

printf("\n"); 

} 

retum(brk); 

} 

int  read_values(void) 

{  . 

int  brk=0,iline, chuff, i_occur,jjp,vaIues_number; 
char  ibuff[256]; 

FILE  *flp  vaiues; 

/***  Reads  in  parameter  values  from  a  file***/ 
printf(" input  file  from  which  to  read  parameter  values:  '*); 
gets{values_filename); 
if  (values_filename[0]  !=  \0') 

{ 

do 

{ 

printf(”input  run  number  of  desired  values  (values_number>0):  "); 
sscanf(gets(ibuff),’'%d”,&values_number); 

}  while  (values_number<  1 ); 

if((fp_values  =  fopen(vaIues_filename;'r”)  )!=  NULL) 

printf(”Searching  file  %s  for  %d  occurrence  of  B10[l]  stiing.\n”, 
values_filename,values_number); 

^^KLoccur=  1  ;i_occur<==values_number;i_occur-H-) 

for(iline-l  ;((fgets(ibufr;  1 00,fp_values)!=NULL)&&!(sscanf(ibuff,"B10[  ll%c” 

&cbuff)));  il  ine++ ) ; 
tf(!{sscanf(ibuff;"B10[  1  ]%c",&cbuff))) 

if(i_occui^l) 

printflj  The  string  B10[l]  was  not  found  in  file  %s  \n”,values_filename); 
else  ~ 

printf("Only  %d  occurrences  of  B10[l]  string  were  found\n”,i_occur); 
brk=l; 

} 

else 

{ 

printfC’Reading  file  %s  starting  at  line  %d\n",values_filename,iiine+l)- 
fgets(ibuff255,fpvalues); 

sscanf(ibuff,"%f  %f  %f  %lf  %lf  %d  %d  %lf  %lf  %lf\n’’,&BlochO[l],&BlochO[2], 
&Bloch0[3],&phiD,&eps,&arun,&Ndetune,&detunei,&detuneN,&inhomo’HW) 

NN=Ndetune-l; 
detunel  *=PI2; 
detuneN  *=  PI2; 
inhomo  HW  *=  PI2; 
arun==0; 

fgets(ibufF,255,:Q3_values); 
fgets(ibufF,255,fj3_values); 
sscanf(ibuff,”%d'\&Npulse); 
if  ((Npulse  >  maxNpuIse)|l(Npuise<l)) 

printf(’'Npulse  =  %d  is  <1  or  >maxNpuise=  %d  \n", Npulse, maxNpulse); 
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brk=2; 

} 

for(jjp=l;((jjp<=Npulse)&&(brk=0))jjp-i“+) 

{ 

for(iline=I;((fgets(ibuff,255,fp_values)!=NlJLL)&&!(sscanf(ibufF,"puise#%c'’, 
&cbuff)));iline-H-); 
if  ((sscanf(ibuff/’pulse#%c",&cbufi))) 

{ 

fgets(ibuf52 5 5 ,fp_values); 

sscanf(ibuff,"%d  %d  %d  %ur’,&jp,&caIc_mode[jjp],&data_flag[ijp], 
&data_value[ijp]); 
if  (jp=ijp) 

{ 

fgets(ibufF,255,fp_values), 

fgets(ibuff,255,fjp_values); 

sscanf(ibuff;rolf%lf%lf%lf%lf%lf%lf%lf",&omegaROp],&tstartOp], 

&tstop[jp],&area,&detuneL[jp],&phiO[jp],&chirpBW[jp],&chirp); 

omegaR[jp]*=PI2; 

detuneL[jp]*=PI2; 

phi01jp]*-PI; 

chirpBW0p]*=PI2; 

area*=PI; 

chirp*=PI2; 


else 

{ 

printf("ErTor;jp=%d  does  not  equal  ijp=%d.\n’' jp,j[jp); 
brk=2; 

} 

} 

else 

{ 

printf(”Error:  string  pulse#  not  found  for  jijp=%d.\n'*,j[jp); 
brk=2; 


} 

fgets(ibuf5255,fp_values); 

if((!(sscanf(ibuff;’idet%c^&cbuflO))&&(b^k^ 

{  . 

printflj "Warning:  The  next  line  doesn't  contain  the  string  idetune  \n")' 

} 

} 

fclose(fpvalues); 

} 

else 

{ 

printf("VaIues_file  %s  can  not  be  opened\n",values_filename); 
brk=l; 

} 

switch(brk) 

{ 

case  0: 

printf("Read  successflilAn"); 
break; 
case  1: 

printf("Read  unsuccessflilAn"); 
break; 
case  2: 
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printf("Read  unsuccessful.  Values  were  altered.  Reinitializing  values:\n"); 

save_flag  =  0; 

initvaiuesO; 

break; 


} 

jp=i; 

retum(brk); 


OBE_BS.C 

#defme  IMAX  1 1 
#define  NUSE  7 
#define  SHRINK  0.95 
#defineGROW  1.2 
#define  NVAR  3 
#include  <math.h> 

#include  <stdio.h> 

#include  "bloch.h" 

#include  "nrutil.h" 

double  d[NVAR+l][NUSE+l],x[IMAX+l];  /*  defining  declaration  */ 
void  bsstep(double  y[],  double  dydx[],  double  *x,  double  htry,  double  eps,  double  yscal[], 
double  *hdid,  double  *hnext,  void  (*derivs  )(),  double  A,  double  B,  double  C, 
double  gamma,  double  omegaR) 

int  ij; 

double  xsav,xest,h,errmax,temp, 

double  ysav[NVAR+ 1  ],dysav[NVAR+ 1  ],yseq[NVAR+l  ],yerr[NVAR+ 1]; 

static  int  nseq[IMAX+l]={0,2, 4,6, 8,12, 16,24,32,48, 64,96}; 

void  mmid(),rzextr(),nrerroit); 

h=htry; 

xsav=(*x); 

for  (i=l;i<=NVAR;i-H-)  { 
ysav[i]=y[i]; 
dysav[i]=dydx[i]; 

!***  Take  square  root  of  yscal  (WRB  1/10/92)  ***! 
yscal[i]  =  exp(0.5*log(yscal[i])); 

} 

for  (;;)  { 

for  (i=l;i<=IMAX;i++)  { 

mmid(ysav,dysav,xsav,h,nseq[i],yseq,derivs,A,B,C,gamma,omegaR); 
xest={temp=h/nseq[i],temp*temp); 
rzextr(i,xest,yseq,y,yerr); 
errmax=0.0; 
for  (j=l;j<=NVARj-H-) 
if  (errmax  <  fabs(yerr|]]/yscal[j])) 
errmax=fabs(yerr[j]/yscal[j]); 
errmax  /=  eps; 
if  (errmax  <  1 .0)  { 

*x  +=  h; 

*hdid=h; 

*hnext  =  i=NUSE?  h*  SHRINK  :  i==NUSE-l? 

h*GROW  ;  {h*nseq[NUSE-l])/nseq[i]; 
return; 

) 

} 

h  *=  0.25; 

for  (i=l;i<=(IMAX-NUSE)/2;i++)  h  /=  2.0; 
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if  ((*x+h)  =  (*x))  nrerrorC’Step  size  underflow  in  BSSTEP"); 


void  inmid(double  y[],  double  dydx[],  double  xs,  double  htot,  int  nstep,  double  yout[], 
^  void  (*derivs  )(),  double  A,  double  B,  double  C,  double  gamma,  double  omegaR) 

int  n,i; 

double  x,swap,h2.h,ym[NVAR+l],yn[NVAR+l], 
h=htot/nstep; 
for  (i=l;i<=NVAR:i-H-)  { 
ym[i]=y[i]; 
yn[i]=y[i]+h*dydx[i], 

} 

x=xs+h; 

(*derivs)(x,yn,yout,A,B,C,gamma,omegaR)- 

h2=2.0*h; 

for  (n=2;n<==nstep;n-K-)  { 
for  (i=l;i<=NVAR;i++)  { 
swap=ym[i]+h2  *yout  [i] ; 
ym[i]=yn[i]; 
yn[i]=swap; 

} 

X  +-  h; 

(*derivs)(x,yn,yout,A,B,C,gamma,omegaR); 

for  (i=l;i<=NVAR;i-H+) 
yout[i]=0.5*(ym[i]+yn[i]+h*yout[i]); 


OBE  FFT.C 

/* 

”fourl()  replaces  data[]  by  its  discrete  Fourier  transform,  if  isign=l;  or 
replaces  data[]  by  nn  times  its  inverse  discrete  Fourier  transform,  if 
isign=-L  data[]  is  a  complex  array  of  length  nn,  input  as  a  real  array 
data[1...2*nn].  nn  MUST  be  an  integer  power  of  2.  (This  is  not  checked!) 

Reference:  Numerical  Recipes,  pp.  407-4 1 2 
♦/ 

^define  PI2  6.28318530717958 

#define  TINY  1 .  Oe- 1 8  /*  *  *  definition  of  a  tiny  value  *  *  */ 

#defineSQR(a)((a)*(a)) 

#define  SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr 
#define  PI  3.14159265358979323846 
#include  <math.h> 

#include  <stdio.h> 

#include  "bloch  h” 

#include  "random,  h" 
static  long  idum2=(-lL); 
static  int  iduml  — 1; 
static  int  idumg  =  -13; 

int  FFT_output(f!oat  _huge  Bloch_matrix[][NVAR+l],  float  CTsignal[],  double  detunel,  double  detuneN 
^  double  inhomo_HW,  double  noise,  double  phiD,  int  NN,  char  CT_fiIename[], double  gamma) 

int  isign,ift,brk,idetune,NNhalf,NN2; 
double  delta_t,  line  shape,  ddetune,  sigmax,  sigmaxt 
FILE  *fp_CT, 
if  (noise  <  0.0) 
idumg  =  (-13), 

printf("Perfonning  FFT  on  output  and  storing  it  in  %s  .\n",CT_filename); 
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delta_t=PI2/(detuneN-detune  1 ); 

/***  convert  density  matrix  to  output  electric  field  ***/ 

NNhalf=NN  »  1 ; 

NN2=NN  «  1; 

ddetune=(detuneN-detune  1  )/((double)NN); 
for{idetune=NNhalf+ 1  ,ift=  1  ;idetune<=NN+ 1  ,idetune-H-,ift  +=2) 

lineshape=1.0; 
if  (inhomo_HW>TINY) 

lineshape=exp(-SQR((detune  1  +ddetune*(double)(idetune- 1  ))/inhomo  HW)) 
if  (inhomo_HW<(-TINY)) 

lineshape=l  .0/(  1 .0+SQR((detunel+ddetune*(double)(idetune-  l))/inhomo_HW)); 
if  (fabs(lineshape)<TINY)  lineshape  =  0.0; 
lineshape  *=  1  +  noise*ddetune*gasdev(&idumg); 
if  (fabs(lineshape)<TINY)  lineshape  =  0.0; 

CTsignal[ift]  =  (float)lineshape*Bloch_matrix[idetune][2]*Bloch_matrix[idetune][0]; 
CTsignal[ift+I]  =  (float)lineshape*Bloch_matrix[idetune][3]*Bloch_matrix[idetune][0]; 

for(idetune=2,ift=NN+3;ift<=NN2- 1  ;idetune++,ift  +=2) 

lineshape=1.0; 
if  (inhomo_HW>TIN  Y) 

Iineshape=exp(-SQR((detune  1  +ddetune*(double)(idetune- 1  ))/inhomo  HW)) 
if(inhomo_HW<(-TINY)) 

lineshape=l  .0/(  1 .0+SQR((detunel  +ddetune*(double)(idetune- 1  ))/mhomo_HW)); 
if  (fabs(lineshape)<TrNY)  lineshape  =  0.0; 
lineshape  *=  1  +  noise*ddetune*gasdev(&idumg); 
if  (fabs(lineshape)<TINY)  lineshape  =  0.0; 

CTsignal[ift]  =  (float)lineshape*Bloch_matrix[idetune][2]*Bloch_nfiatrix[idetune][0]; 

^  CTsignal[ift+l]  =  (float)lineshape*Bloch_matrix[idetune][3]*Bloch_matrix[idetune][’o]; 

isign=l; 

FFT(C  Tsignal,  NN,  isign ) ; 

/***  Store  CTsignal  out  for  graphit  to  graph  ***/ 
if((  fp_CT  =  fopen(CT_filename,"w")  )=  NULL) 

{ 

printf("CT_File  %s  can  not  be  opened", CT  filename); 
brk=l; 

} 

else 

{ 

sigmax  =0.0; 
sigmaxt  =  0.0; 

for(ift=l;ift<=NN2-l;ift  +=2)  { 

^rintfi:fp_CT,  "%g\t%.5g\t%.  5g  \n", 

(double)(ift- 1  )*0.5*delta_t,(SQR(CTsignal[ift])+SQR(CTsignal[ift+ 1  ]))* 
exp(-(((double)(ift-l)*0.5*delta_t-190)*gamma*2.0)), 
(cos(phiD)*CTsignal[ift+l]-sin(phiD)*CTsignal[ift])* 
exp(-(((double)(ift- 1  )*0. 5*delta_t- 1 90)*gamma)) ); 
if((double)(ift-l)*0.5*delta_t  >  70.0) 
if  ((SQR(CTsignal[ift])+SQR(CTsignal[ift+l]))>  sigmax)  { 
sigmax  =(SQR(CTsignal[ift])+SQR(CTsignal[ift+l])); 
sigmaxt  =(double)(ift- 1  )*0. 5  *delta  t; 

} 


fclose(fp_CT). 

printf("sigmaxt  =  %12.5g  simmax  =  %12.5g  \n",sigmaxt,sigmax); 
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retum(brk); 


} 

OBE_INT.C 

/*This  program  is  to  be  used  in  conjunction  with  bloch.c  for 
integrating  the  Optical  Bloch  Equations. 

It  is  based  on  odeint.c,  except  that  intemediate  value  of  the  integration 
are  no  longer  stored  in  xp  and  yp.*/ 

#defineMAXSTP  10000 
#define  TINY  l  .Oe-30 
#defineNVAR  3 
?Mefine  PI  3.14159265358979 
#define  PI2  6.28318530717958 
#include  <stdlib.h> 

^include  <bios.h> 

#include  <conio.h> 

#include  <io.h> 

^include  <stdio.h> 

^include  <math.h> 

#include  "complex. h” 

^include  "bloch.h" 

#include  "nrutil.h" 

#include  "stdio_x.h" 
static  int  idumg2=-13; 

int  integrate_pulse(float  _huge  BIoch_matrix[][NVAR+l],  FILE  *fp_save,  int  jpulse,  int  caic  mode, 
double  eps,  double  omegaR,  double  tstart,  double  tstop,  double  detuneL,  double  phiO, 
double  chirpBW,  double  gamma,  double  detune!,  double  ddetune,  int  Ndetune, 
int  printout_flag,  int  save  flag) 

int  nbad,nok,idetune,ivar,brk=0; 
unsigned  key  =  0; 

double  hl,hinin,Bloch[NVAR+l],deltaA,deltaB,A,B,C,detune,detuneAL,tau,area,chirp; 

tau=tstop-tstart; 
area=  omegaR*tau; 

chirp=(fabs(chirpBW)>TINY)?(chirpBW/tau):0.0; 
printf("%6s%4s\n","pulse#",''calc"); 
printf("%6d  %4d\n",jpulse,calc_mode); 
if  (calc_mode==5) 
printf("gamma  =  %lf\n", gamma); 
printf(’’%12s  %8s  %8s  %10s  %7s  %5s  %7s  %10s  \n", 

"omegaR","tstart",  "tstop", "area",  "detuneL",'’phiO'',"chirpBW",”chirp") 
printf("%12.6g  %8,3f%8.3f%10.3g  %7.3f%5.3f  %7.3f  %10,3g\n", 
(omegaR/PI2),tstart,tstop,(areayPI),(detuneL/PI2), 

(phiO/PI),(chirpBW/PI2),(chirp/PI2)); 
if  (printout  flag  =  1) 

printf("%4s  %7s  %9s  %10s  %10s  %9s  %4s  %4s\n","idet", 

"detune","Bloch[l]","Bloch[2]","Bloch[3]","Bloch[0]","nok'’,"nbad"); 

if  (save  flag  ==  1 ) 

{ 

fprintf(fp_save,"%6s\t%4s\t%4s\t%10s\n","pulse#","calc","data","data_value")- 
fprintf(^_save,"%6d\t%4d\t%4d\t%10d\n",jpulse,calc_mode,0,0); 
fprintf(f^_save,"%  1 2s\t%8s\t%8s\t%  1 0s\t%7s\t%5s\t%7s\t%  1 6s 
"omegaR","tstart'',  "tstop",  "area","detuneL","phiO","chirpBW", "chirp"); 
fprintf(fp_save,"%12.6g\t%8.3f\t%8.3fvt%I0.3g\t%7.3flt%5.3f\t%7.3f\t%io,3g\n" 
(omegaR/PI2),tstart,tstop,(area/PI),(detuneL/PI2), 

(phiO/PI),(chirpBW/PI2),(chirp/PI2)); 

fprintf(fp_save,"%4s\t%7s\t%9s\t%10s\t%10s\t%9s\t%4s\t%4s\n", 

"idet","detune","Bloch[l]","Bloch[2]","Bloch[3]","Bloch[0]","nok","nbad"); 
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} 

cietune=detunel; 
detxineAL=  detune  1  -  detuneL; 

B=  -detuneAL; 

C=0.5*chirp; 

A=  phiO  -  detuneAL  *tstart; 
deitaA  =  -ddetune*tstart: 
deltas  =  -ddetune; 
switch(calc_mode) 

{ 

case-1;  /***  do  nothing  ***/ 
break; 

case  0:  /***  Runge-Kutta  integration  ***/ 

hmin=0.0; 
hl-0.1; 
break; 

case  1:  /***  Bulirsch-Stoer  integration  with  constant  amplitude***/ 
hmin=0.0; 
hl=2.0, 
break; 

case  2:  /***  Analytic  solution  ***/ 
nok  =  1 ; 
nbad  0; 
break; 

case  3 :  /*  *  *  coherence  loss  *  *  */ 
break; 

case  4;  /***  gating  pulse  ***/ 
break; 

case  5:  /***  Bulirsch-Stoer  integration  with  decay***/ 

hmin-0.0; 
hl-2.0; 
break; 

case  6:  /***  Bulirsch-Stoer  integration  with  shaped  pulses***/ 

hmin=0.0; 
hl=2.0; 
break; 

case  7;  /***  gating  pulse  and  noise***/ 
break; 

} 

for(idetune=l;((idetune  <=  Ndetune)&&(brk=0));  idetune++) 

/***  set  Bloch  to  old  value  of  Bloch  matrix  ***/ 

for(ivar=0,ivar  <—  NVAR;  ivar“H-)  Bloch[ivar]=(double)  Bloch_matrix[idetune][ivar]; 
switch(calc_mode) 

{ 

case-1:  /***  do  nothing  ***/ 
break; 

case  0:  /***  Runge-Kutta  integration  ***/ 

if  (omegaR>TrNY) 

OBEint(Bloch,0.0,tau,eps,hl,hmin,&nok,&nbad,nodecay,A,B,C,0.0, 

omegaR^rkqc); 

break; 

case  1 :  /***  Bulirsch-Stoer  integration  with  constant  amp***/ 

if  (omegaR>TrivJY) 

OBEint(BIoch,0.0,tau,eps,hl,hmin,&nok,&nbad,nodecay,A,B,C,gamma, 

omegaR^bsstep),  /*  for  constant  omegaRt  gamma<=>  0.0  moe  12/19/96  */ 
break; 

case  2;  /***  Analytic  solution  ***/ 


33 


if  (omegaR>TINY) 

SquarePulse(Bloch,tau,A,-B,C,omegaR); 

break; 

case  3:  /***  coherence  loss  ***/ 

Bloch[2]=0,0; 

Bloch[3]=0.0; 

break; 

case  4:  /***  gating  pulse  ***/ 

Bloch[2]=0.0; 

Bloch[3]=0.0; 

Bloch[0]  ♦=  Bloch[l]; 

Bloch[l]=l,0; 

break; 

case  5:  /*♦*  Bulirsch-Stoer  integration  with  decay***/ 

if  (omegaR>TINY) 

OBEint(Bloch,0,0,tau,eps,hl,hniin,&nok,&nbad,decay,A,B,C,gamma, 

oniegaR,bsstep);  /*assumes  constant  amplitude  pulse*/ 
break; 

case  6:  /***  Bulirsch-Stoer  integration  with  shaped  pulse***/ 

if  (oniegaR>TINY) 

OBEint(Bloch,0.0,tau,eps,hl,hniin,&nok,&nbad,nodecay,A,B,C,(PI2/tau), 

omegaR,bsstep);  /*for  shaped  pulses*/ 
break; 

case?:  /***  gating  pulse  and  noise***/ 

Bloch[2]=0.0; 

Bloch[3]=0.0; 

Bloch[0]  *=  Bloch[l]; 

Bloch[0]  +=  omegaR*ddetune*gasdev(&idumg2); 

Bloch[l]=l,0; 

break; 

} 

/***  store  new  value  in  Bloch  matrix  ***/ 
for(ivar=0;ivar  <=  NVAR;  ivar++) 

{ 

modf(  1 00000000. 0*Bloch[ivar],&Bloch[ivar]);  **  *  1  e-8  accuracy*  *  * 
Bloch_matrix[idetune][ivar]  =  Bloch[ivar]*0.00000001;  */ 

^  Bloch_matrix[idetune][ivar]  =  (float)Bloch[ivar]; 

if  (printoutflag  =  1 ) 

printf("%4d%7.4f%9.7f%10.7f%10.7f%9,7f%4d%4d\n",idetune,(detune/PI2), 
Bloch_matrix[idetune][l  ],Bloch_matrix[idetune][2], 
Bloch_matrix[idetune][3],Bloch_matrix[idetune][0]’nok,nbad); 
if((save_flag=  1  )&((printout_flag  !=  2)||(idetune=l))) 
t)rintfl:fp_save,"%d\t%.4g\t%g\t%g\t%g\t%g\t%d\t%d\n", 
idetune,(detune/PI2),BIoch_matrix[idetune][  1  ],Bloch_matrix[idetune]  [2], 
Bloch_matrix[idetune][3],Bloch_matrix[idetune][0],nok,nbad); 

A  +=  deltaA; 

B  +=  deltaB; 
detune  +=  ddetune; 

if  C.bios_keybrd(_NKEYBRD_READY)) 

key  =  _bios_keybrd(  NKEYBRDJREAD  ); 
if  ((key  &  OxOOff)  =  'b') 

{ 

brk  =  1 ; 

printf("break\n"); 
if  (save_flag  =  1 ) 
fprintf(fp_save,"%6d\n",- 1 ); 


34 


} 

} 

}  /***  end  for(ident=  ..)  ***/ 
retum(brk); 

} 

void  OBEint(double  ystart[],  double  xl,  double  x2,  double  eps,  double  hi,  double  hmin, 
int  *nok,  int  *nbad,void  (*derivs  )(),  double  A,  double  B, 

^  double  C,  double  gamma,  double  omegaR,  void  (*  stepper  )()) 

int  nstp,i; 

double  x,hnext,hdid,h; 

double  yscal[NVAR+ 1  ],y[NVAR+ 1  ],dydx[NVAR+ 1  ]; 

void  nrerrorO; 

x=xl; 

h={x2  >  xl)  ?  fabs(hl)  :  -fabs(hl); 

*nok  =  (*nbad)  =  0; 

for  (i=l;i<=NVAR;i-i-+)  y[i]=ystart[i]; 

for  (nstp=l;nstp<=MAXSTP;nstp++)  { 

(*derivs)(x,y,dydx,A,B,C, gamma,  omegaR); 
if  ((x+h-x2)*(x+h“xl)  >  0.0)  h=x2-x; 
for  (i=l;i<=NVAR;i++) 

/***  Use  absolute  accuracy  scaled  to  step  size  (WRB  1/10/92)  ***/ 
yscal[i]=fabs(h/(x2-x  1 )); 

(*stepper)(y,dydx,&x,h,eps,yscal,&hdid,&hnext,derivs,A,B,C,gamma,omegaR); 

if  (hdid  =“  h)  -H-(*nok);  else  ^(*nbad); 
if  ((x-x2)*(x2-xl)  >=  0.0)  { 
for  (i=l;i<=NVAR;i-H-)  ystart[i]=y[i], 
return; 

) 

if  (fabs(hnext)  <=  hmin)  nrerror("Step  size  too  small  in  OBEINT”); 
h=hnext; 

} 

nrerror(”Too  many  steps  in  routine  OBEINT”); 
return; 

} 

void  nodecay(double  t,  double  BIoch[],  double  dBlochdt[],  double  A,  double  B,  double  C, 
double  gamma,  double  omegaR) 

double  phit,omegaRt,al,a2,a3; 
if{gamma=0,0)  omegaRt  =  omegaR; 

else  omegaRt  =  omegaR*(  1 .0-cos(PI2*t*gamma)*cos(PI2*t*gamma))'  /*PI2  et  al 
moe  12/19/96*/ 
phit=  A  +  B*t  +  C*t*t; 
al=  omegaRt*sin(phit); 
a2=  omegaRt*cos(phit); 
a3=  0.5 -Bloch[l]; 

dBlochdt[l]=  (a2*Bloch[3]  -  al*Bloch[2]); 
dBlochdt[2]=  -al*a3; 
dBlochdt[3]=  a2*a3 

}  . 

void  decay(double  t,  double  Bloch[],  double  dBlochdt[],  double  A,  double  B,  double  C, 
double  gamma,  double  omegaR) 

double  phit,omegaRt,al,a2,a3; 
phit=  A  +  B*t  +  C*t*t; 
al=  omegaR*sin(phit); 
a2=  omegaR*cos(phit); 
a3=  0.5-Bloch[l]; 
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dBlochdt[l]=  (a2*Bloch[3]  -  al*Bloch[2]); 
dBlochdt[2]=  -al  *a3-gamma*Bloch[2], 
dBlochdt[3]=  a2*a3-gamma*Bloch[3]; 

OBE_RK.C 

^define  PGROW  -0.20 
#define  PSHRNK  -0  25 

#define  FCOR  0.06666666  /*  1  0/15  0  */ 

#define  SAFETY  0.9 
#define  ERRCON  6.0e-4 
#defineNVAR  3 
#include  <math.h> 

#include  <stdio.h> 

#include  "nrutil.h" 

#include  "bloch.h" 

void  rkqc(double  y[],  double  dydx[],  double  *x,  double  htry,  double  eps,  double  yscal[], 
double  *hdid,  double  *hnext,  void  (*derivs  )(), 
double  A,  double  B,  double  C,  double  gamma,  double  omegaR) 

int  i; 

double  xsav,hh,h,temp,errmax; 

double  dysav[NVAR+ 1  ],ysav[NVAR+ 1  ],ytemp[NVAR+l  ]; 

void  rlc4(),nreiTor(); 

xsav=(*x); 

for  (i=l;i<=NVAR;i++)  { 
ysav[i]=y(i]; 
dysav[i]=ciydx[i]; 

/***  Take  square  root  of  yscal  (WRB  1/10/92)  ***/ 
yscal[i]  =  sqrt(yscal[i])*16.0; 

} 

h=htry; 
for  (;;)  { 
hh=0.5*h; 

rk4(ysav,dysav,xsav,hh,ytemp,derivs,A,B,C,gamma,omegaR); 

*x=xsav+hh; 

(*derivs)(*x,ytemp,dydx,A,B,C,gamma,omegaR); 
rk4(ytemp,dydx,*x,hh,y,derivs,A,B,C, gamma,  omegaR); 

*x=xsav+h; 

if  (*x  =  xsav)  nrerrorC'Step  size  too  small  in  routine  RKQC"); 
rk4(ysav,dysav,  xsav,  h,ytemp,derivs,  A,B,  C,gamma,  omegaR); 
eiTmax=0.0; 

for  (i=l;i<=NVAR;i-H-)  { 
ytemp[i]=y[i]-ytemp[i]; 
temp=fabs(ytemp[i]/yscal[i]); 
if  (errmax  <  temp)  ernnax=temp; 

} 

errmax  /=  eps; 
if  (errmax  <=  1 .0)  { 

*hdid=h; 

*hnext=(ernnax  >  ERRCON  ? 

SAFETY*h*exp(PGROW*log(errmax))  :  4.0*h); 
break; 

} 

h=SAFETY*h*exp(PSHRNK*log(errmax)); 
for  (i=l;i<=NVAR;i-H-)  y[i]  -i-=  ytemp[i]*FCOR; 
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void  rk4(double  y[],  double  dydx[],  double  x,  double  h,  double  yout[], 
void  (*derivs  )(),  double  A,  double  B,  double  C,  double  gamma, 
double  omegaR) 

{ 

int  i; 

double  xh,hh,h6,dym[NV AR+ 1  ],dyt[NVAR+ 1  ],yt [NV AR+ 1  ] • 

hh=h*0.5; 

h6-h/6.0; 

xh=x+hh; 

for  (i=l;i<=NVAR;i++)  yt[i]=y[i]+hh*dydx[i]; 
(*derivs)(xh,yt,dyt,A,B,C,gamma,omegaR); 
for  (i=l;i<=IWAR;i-H“)  yt[i]=y[i]+hh*dyt[i]; 
(*derivs)(xh,yt,dym,A,B,C,gamma,  omegaR); 
for  (i-l;i<=NyAR;i+-H)  { 
yt[i]=^[i]+h*dym[i]; 
dym[i]  +=  dyt[i], 

} 

(*derivs)(x+h,yt,dyt,A,B,C,gamma,omegaR); 
for  (i=l;i<=NVAR;i++) 
yout[i]=7[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]); 


OBE^SQX 

#include  <stdio.h> 

#include  <math.h> 

^include  "complex. h" 

#include  "bloch.h" 

#include  "nrutil.h" 

#defineSQR(a)((a)*(a)) 

void  SquarePulse(double  y[],  double  tau,  double  phi,  double  delta,  double  C,  double  rabi) 
double  chi,theta,theta_2; 

fcomplex  rho01_start,rho01_final,rho00_start,rho00_final,Cphase,Ctempl,Ctemp2,Ctemp3; 

void  nrerrorO; 

if(CI=0.0) 

nrerror("Chirp  does  not  equal  zero"); 
else 
{ 

rho00_start  =  Complex(y[l],0.0); 
rho01_start  =  Complex(y[2],y[3]), 
chi  =  sqrt(SQR(delta)  +  SQR(rabi)); 
theta  =  chi*tau; 
theta_2  =  0.5*theta; 

Cphase  =  Complex(0,phi); 

Ctempl=Cmul( 

Compiex((SQR(cos(theta_2))-SQR(delta*sin(theta_2)/chi)), 

(delta*  sin(thetaychi)), 
rho01_start); 

Ctemp2=Cmul( 

RCmul( 

SQR(rabi  *  sin(theta_2)/chi), 

Cexp(RCmul(2 . 0,  Cphase))), 

Conjg(rho01  start) ); 

Ctemp3=Cmul( 

Complex( 

delta*rabi*SQR(sin(theta_2)/chi),  - 
rabi*0.5*sin(theta)/chi) , 

Cexp(Cphase)  ); 
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rhoO  1  _final=C  mul( 

Cadd( 

Cadd( 

Ctempl, 

Ctenip2  ), 

RCmul( 

(2.0*rho00_start.r-  1.0), 

Ctemp3)), 

Cexp(Complex(0.0,-delta*tau) ) ); 

rho00_finaI.r=rho00_start.r  +  (1 .0-2.0*rho00_start.r)*SQR(rabi*sin(theta_2)/chi) 
+(rabi*sin(theta)/chi)*(cos(phi)*rho01_start.i  -  sin(phi)*rho01_start.r) 
+2.0*rabi*delta*SQR(sin(theta_2ychi)* 

(sin(phi)*rho01_start.i  +  cos(phi)*rho01_start.r); 
y[2]=rho01_final.r; 
y[3]=rho01  final,  i; 
y[l]=rho00_final.r; 

} 

return; 

} 
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